// Programa que calcula una norma para cada una de las     |
// filas de una matriz bidimensional representada mediante |
// un formato de almacenamiento comprimido.                |
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <pmmintrin.h>

//TODO Mejora de que no falle cuando el numero de filas no es multiplo de 4


#define INDEX(i,j) [(i)*ncols+j]

void inicializacion(float **pvalores, float **pnorma, int nfilas, int ncols, int nfilasaux, int ncolsaux);
__m128 hadd(__m128 a, __m128 b);


int main( int argc, char *argv[] ) {
	
	float *valores; 
	float *norma;
	float normafila;
	int nfilas, nfilasaux ,ncols, ncolsaux;
	int i, j;
	__m128 reg_valores, reg_valores2, reg_sum, reg_aux1, reg_aux2;
	struct timeval t0, t1, t;
	
	float *aux;//Auxiliar para guardar los cuatro datos multiplicados en cada iter.
	
	
	if (argc != 3) {
		printf ("Uso: norma nfilas ncols \n");
		exit(0);
	}
	nfilas = atoi(argv[1]);
	ncols = atoi(argv[2]);

        nfilasaux = nfilas;
	/* Nfilas ahora es multiplo de cuatro (redondeado hacia arriba)*/        
        if (nfilas%4)
		nfilas = nfilas + (4-(nfilas%4));

	ncolsaux = ncols;	
	/* Ncols ahora es multiplo de cuatro (redondeado hacia arriba)*/
	if(ncols % 4)
		ncols = ncols + (4-(ncols%4));

	inicializacion(&valores, &norma, nfilas, ncols, nfilasaux, ncolsaux);
	assert (gettimeofday (&t0, NULL) == 0);
			
	aux = (float *)_mm_malloc(4 * sizeof(float), 16);
	for (i=0; i<nfilas; i+=4) {
		reg_sum = _mm_setzero_ps();
		for (j=0; j<ncols; j+=4) {		
								
		        reg_valores = _mm_load_ps(&(valores INDEX(i,j)));
		        reg_valores = _mm_mul_ps(reg_valores, reg_valores);
		        
			reg_valores2 = _mm_load_ps(&(valores INDEX(i+1,j))); 
			reg_valores2 = _mm_mul_ps(reg_valores2, reg_valores2);		       
		       
			reg_aux1 = hadd(reg_valores, reg_valores2);
					
						
		        reg_valores = _mm_load_ps(&(valores INDEX(i+2,j)));
		        reg_valores = _mm_mul_ps(reg_valores, reg_valores);
		        
			reg_valores2 = _mm_load_ps(&(valores INDEX(i+3,j))); 
			reg_valores2 = _mm_mul_ps(reg_valores2, reg_valores2);		        

			reg_aux2 = hadd(reg_valores, reg_valores2);			
						
			reg_aux2 = hadd(reg_aux1, reg_aux2);		

			reg_sum = _mm_add_ps(reg_sum, reg_aux2);													
		}

		reg_sum = _mm_sqrt_ps(reg_sum);
		
		_mm_store_ps(aux, reg_sum);
		for(j=0; j<4; j++)
		        if (i+j < nfilasaux)   
			        norma[i+j] = aux[j];
	}	
	
	
	assert (gettimeofday (&t1, NULL) == 0);
	timersub(&t1, &t0, &t);

	for(i=0;i<5;i++) printf("Norma[%d] = %.2f \n",i,norma[i]);
        for(i=nfilasaux-5;i<nfilasaux;i++) printf("Norma[%d] = %.2f \n",i,norma[i]);
	
	printf ("nfilas      = %ld\n", nfilasaux);
	printf ("ncols       = %ld\n", ncolsaux);
	printf ("Tiempo      = %ld:%ld(seg:mseg)\n", t.tv_sec, t.tv_usec/1000);
	
	free(valores); 
	free(norma);
	
	return 0;
}


void inicializacion(float **pvalores, float **pnorma, int nfilas, int ncols, int nfilasaux, int ncolsaux)
{
	int i, j;

	
	if( ((*pvalores)=(float *) _mm_malloc(sizeof(float)*nfilas*ncols,16)) == NULL ) {
		printf ("error: memoria insuficiente para valores\n");
		exit(-1);
	}

	if( ((*pnorma)=(float *) malloc(sizeof(float)*nfilasaux)) == NULL ) {
		printf ("error: memoria insuficiente para norma\n");
		exit(-1);
	}

	for (i=0; i<nfilas; i++) {	
	        if (i<nfilasaux)
		        for (j=0; j<ncols; j++) {
			        if(j<ncolsaux)
				        (*pvalores) INDEX(i,j) = (float) ( (float) random() / (float) RAND_MAX );
			        else
				        (*pvalores) INDEX(i,j) = 0;  // inicializamos a cero las columnas añadidas a mayores
		        }
		else 
                        for (j=0; j<ncols; j++) 
				(*pvalores) INDEX(i,j) = 0;  // inicializamos a cero las filas añadidas a mayores
		        
	}
}

__m128 hadd(__m128 a, __m128 b)
{
	__m128 aux1, aux2;
	
	aux1 = _mm_shuffle_ps(a, b, _MM_SHUFFLE(2,0,2,0));
	aux2 = _mm_shuffle_ps(a, b, _MM_SHUFFLE(3,1,3,1));

	return _mm_add_ps(aux1, aux2);
}

